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Abstract 

In  this  paper  we  consider  a  strongly  coupled  (monolithic)  fluid  structure  interaction  framework  for  incom¬ 
pressible  flow,  as  opposed  to  a  loosely  coupled  (partitioned)  method.  This  requires  solving  a  single  linear 
system  that  combines  the  unknown  velocities  of  the  structure  with  the  unknown  pressures  of  the  fluid.  In 
our  previous  work,  we  were  able  to  obtain  a  symmetric  formulation  of  this  coupled  system;  however,  it  was 
also  indefinite,  making  it  more  difficult  to  solve.  In  fact  in  practice  there  have  been  cases  where  we  have 
been  unable  to  invert  the  system.  In  this  paper  we  take  a  novel  approach  that  consists  of  factoring  the 
damping  matrix  of  deformable  structures  and  show  that  this  can  be  used  to  obtain  a  symmetric  positive 
definite  system,  at  least  to  the  extent  that  the  uncoupled  systems  were  symmetric  positive  definite.  We  use 
a  traditional  MAC  grid  discretization  of  the  fluid  and  a  fully  Lagrangian  discretization  of  the  structures  for 
the  sake  of  exposition,  noting  that  our  procedure  can  be  generalized  to  other  scenarios.  For  the  special  case 
of  rigid  bodies,  where  there  are  no  internal  damping  forces,  we  exactly  recover  the  system  of  [4]. 


1.  Introduction 

Fluid  structure  interaction  has  been  of  increasing  interest  to  the  computational  fluids  community,  due 
in  part  to  the  ability  to  address  more  computationally  expensive  problems  with  improved  hardware.  Re¬ 
searchers  have  found  differing  approaches  useful  for  qualitatively  different  regimes  of  fluid  structure  interac¬ 
tion.  Applications  may  involve  low  or  high  structural  deformation.  Problems  involving  low  deformation  can 
be  approached  with  aligning  fluid  meshes  using  for  instance  Arbitrary  Lagrangian-Eulerian  (ALE)  schemes. 
This  has  been  quite  successful  for  studying  phenomena  such  as  airplane  wing  flutter,  see  e.g.  [10].  High 
structural  deformation  is  sometimes  also  approached  using  ALE  techniques,  but  frequent  remeshing  and  the 
averaging  it  requires  can  lead  to  reduced  accuracy  and  increased  computational  cost.  More  often,  methods 
where  the  fluid  and  structure  meshes  are  not  aligned  are  employed  (e.g.  the  body  of  work  originating  from 
[23]).  In  this  case  the  interaction  between  fluid  and  structure  is  usually  massaged  via  penalty  methods  or 
Lagrange  multiplier  formulations. 

Anecdotal  evidence  suggests  that  problems  where  either  the  fluid  or  the  structure  dominates  the  inter¬ 
action  are  inherently  more  stable  than  those  where  the  effect  is  more  balanced.  For  example,  it  is  more 
straightforward  to  solve  problems  where  a  kinematic  structure  drives  the  fluid  or  a  light  structure  is  pas¬ 
sively  advected  by  fluid,  as  opposed  to  problems  where  structures  and  fluids  have  nearly  equal  densities  (see 
e.g.  [7]).  Partitioned  coupling  approaches  are  typically  explicit  and  staggered,  for  instance  first  imposing 
the  fluid  pressures  as  a  force  on  the  structure  and  then  solving  for  the  fluid  using  the  structure  velocities 
as  boundary  conditions  in  a  second  step.  These  approaches  can  also  be  iterated  within  a  given  time  step 
for  increased  stability,  noting  that  in  the  limit  if  one  converges  one  obtains  a  monolithic  (albeit  expensive) 
approach.  Other  approaches  construct  strongly  coupled  systems  and  then  solve  them  in  one  of  several  ways. 
We  follow  a  monolithically  coupled  approach,  in  which  we  model  the  fluid  structure  interaction  by  formu¬ 
lating  and  solving  a  single  coupled  system  of  equations;  however,  we  do  this  in  a  way  which  is  more  efficient 


*{avir,cas43,fedkiw}@cs. stanford.edu,  Stanford  University 


Preprint  submitted  to  Elsevier 


August  9,  2010 


Report  Documentation  Page 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 


1.  REPORT  DATE 

09  AUG  2010 


2.  REPORT  TYPE 


4.  TITLE  AND  SUBTITLE 

A  symmetric  positive  definite  formulation  for  monolithic  fluid  structure 
interaction 

6.  AUTHOR(S) 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Stanford  University, 353  Serra  Mall  Room  207, Stanford, CA, 94305 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


3.  DATES  COVERED 

00-00-2010  to  00-00-2010 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

5d.  PROIECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

In  this  paper  we  consider  a  strongly  coupled  (monolithic)  fluid  structure  interaction  framework  for  incom¬ 
pressible  flow,  as  opposed  to  a  loosely  coupled  (partitioned)  method.  This  requires  solving  a  single  linear 
system  that  combines  the  unknown  velocities  of  the  structure  with  the  unknown  pressures  of  the  fluid.  In 
our  previous  work,  we  were  able  to  obtain  a  symmetric  formulation  of  this  coupled  system;  however,  it  was 
also  indefinite,  making  it  more  difficult  to  solve.  In  fact  in  practice  there  have  been  cases  where  we  have 
been  unable  to  invert  the  system.  In  this  paper  we  take  a  novel  approach  that  consists  of  factoring  the 
damping  matrix  of  deformable  structures  and  show  that  this  can  be  used  to  obtain  a  symmetric  positive 
definite  system,  at  least  to  the  extent  that  the  uncoupled  systems  were  symmetric  positive  definite.  We  use  a 
traditional  MAC  grid  discretization  of  the  fluid  and  a  fully  Lagrangian  discretization  of  the  structures  for 
the  sake  of  exposition,  noting  that  our  procedure  can  be  generalized  to  other  scenarios.  For  the  special  case 
of  rigid  bodies,  where  there  are  no  internal  damping  forces,  we  exactly  recover  the  system  of  [4]. 

15.  SUBIECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OF 

ABSTRACT 

OF  PAGES 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

25 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


as  well  as  more  likely  to  converge  than  simply  iterating  the  partitioned  approach  to  convergence  in  a  simple 
Gauss-Seidel  manner. 

Our  approach  allows  the  use  of  existing  discretizations  for  both  the  fluid  and  the  structure  and  requires 
only  minimal  invasion  into  existing  code  bases  for  solving  the  system.  However,  we  do  require  that  the 
fluid  and  structure  solvers  have  certain  properties,  albeit  we  will  argue  that  these  are  properties  that  one 
would  desire  any  physically  constituted  system  to  obey.  The  current  formulation  does  assume  that  the  fluid 
and  structure  systems  are  synchronized  in  time  in  a  manner  that  allows  us  to  solve  for  the  coupling  at  the 
synchronization  points,  although  it  is  not  inconceivable  that  one  could  apply  asynchronous  subiteration  to 
each  system. 

The  method  that  we  propose  in  this  paper  is  similar  to  the  previous  work  in  [25]  addressing  an  implicit 
method  for  two-way  coupling  of  incompressible  fluids  and  structures,  with  the  chief  difference  being  that  the 
linear  system  presented  in  [25]  was  indefinite.  Our  main  contribution  is  the  observation  in  Section  5.2  that 
the  structure  damping  operator  can  be  factored  and  used  to  perform  a  change  of  variables  that  allows  us  to 
symmetrize  the  coupled  system  while  maintaining  positive  definiteness.  We  examine  the  properties  of  the 
resulting  linear  system  and  demonstrate  convergence  of  the  method  on  several  example  problems  from  the 
literature. 


2.  Fluid  Equations 

We  address  the  Navier-Stokes  equations  for  incompressible  flow, 

p(j^  +  (u '  v)u)  =  + f  w 

and 

V  •  u  =  0,  (2) 

in  the  presence  of  structures.  Here,  u  is  the  fluid  velocity,  p  is  the  pressure,  p  is  the  dynamic  viscosity,  p  is 
the  density,  and  f  is  any  external  body  force  such  as  gravity. 

We  use  a  standard  Marker  and  Cell  (MAC,  [16])  grid  discretization  for  the  fluid.  When  solving  in 
the  presence  of  structures,  some  faces  may  not  contain  valid  velocity  samples  because  they  are  covered  by 
structures.  The  vector  of  fluid  velocity  samples  omitting  any  covered  by  the  structure  is  denoted  by  u. 
Pressure  samples  are  defined  at  cell  centers.  The  divergence  operator  is  defined  to  act  on  velocities  at  cell 
faces  to  yield  divergences  at  cell  centers  using  the  standard  six-point  stencil  in  three  spatial  dimensions,  and 
the  gradient  operator  is  defined  as  the  negative  transpose  of  the  divergence. 

We  apply  the  projection  method  of  [8]  to  solve  the  discretized  form  of  Eqs.  (1)  and  (2),  making  a  first  order 
approximation  in  time.  This  splits  the  time  integration  into  the  computation  of  an  intermediate  velocity 
ignoring  the  pressure  terms 

u*  =  un  -  At(un  ■  V)u"  +  ^V2u  +  — f  (3) 

P  P 

followed  by  enforcement  of  incompressibility  via  an  implicit  pressure  projection 

un+1  =  u*  -  — Vp.  (4) 

P 

We  have  omitted  a  specification  of  time  from  the  viscosity  term,  since  the  viscosity  might  be  applied  explicitly 
or  implicitly.  We  postpone  further  consideration  of  viscosity  until  Section  7,  at  which  point  we  will  include 
it  in  the  coupling  framework.  We  use  a  semi-Lagrangian  advection  scheme  (see  e.g.  [26])  to  compute  the 
convective  terms  in  Eq.  (3).  In  the  presence  of  structures,  this  requires  ghost  cells  or  clipping  since  the 
back-cast  rays  may  intersect  the  solid;  we  use  a  method  based  on  [14]  for  semi-Lagrangian  advection  near 


2 


solid  surfaces.  We  advect  using  the  divergence  free  time  n  velocities.  We  enforce  incompressibility  as  usual 
by  plugging  Eq.  (4)  into  the  discrete  form  of  Eq.  (2)  to  obtain 


V  • 


=  V  •  u 


★ 


which  we  solve  for  p  using  the  conjugate  gradient  method. 


3.  Structure  Equations 

We  utilize  a  traditional  Lagrangian  approach  for  structures.  For  deformable  structures,  we  solve  the 
equations  x  =  v  and  Mv  =  F.  We  take  a  lumped  mass  approach  where  the  quantities  x  and  v  are  the 
vectors  containing  the  position  and  velocity  degrees  of  freedom  for  all  n  of  the  particles.  The  mass  matrix 
M  is  a  3n  x  3n  diagonal  matrix  with  each  3x3  diagonal  block  of  the  form  rrijl,  where  to*  is  the  mass  of 
particle  i.  The  forces  F  =  Fe  +  contain  both  elastic  Fe  (velocity  independent)  and  damping  F^  (velocity 
dependent)  forces.  Purely  elastic  forces  (e.g.,  no  plasticity)  can  be  expressed  in  terms  of  an  energy  potential 
as  Fe  =  —  where  the  energy  potential  </>  depends  only  on  x.  Damping  forces  are  often  chosen  according 
to  a  simple  model  F,j  =  Dv  that  is  linear  in  velocities  and  has  symmetric  negative  semidefinite  D.  This  is 
true,  for  example,  for  the  damping  applied  to  a  simple  spring  and  Rayleigh  damping. 

One  particularly  simple  and  stable  way  of  evolving  these  equations  is  via  backward  Euler.  We  wish  to 
solve 

(  x"  +  A  tvn+1  \ 
y vn  +  AtM_1Frt+1y  ' 

Taking  the  first  order  approximation  F,i+1  =  F"  +  ^r(xn+1  —  x")  +  ^-(v"+1  —  v”)  this  becomes 

x"  +  Atv"+1 

v11  +  AtM_1(Fn  +  ^(xn+1  -  x")  +  ^(vn+1  -  v”)) 

Using  the  first  equation  to  eliminate  xn+1  in  the  second  equation  yields 

fy pn  f)T?n 

vn+1  =  v"  +  AtM_1(F”  +  — —  Afv”+1  +  — —  (v"+1  -  v”)) 

ax  av 

(Dun  BV'n  \  r)V'n 

I  -  Ativr1— - Ai2M_1  — —  v”+1  =  vn  -  AfM-1  — — v"  +  AtM-1Fn  (5) 

av  ax  )  av 

From  this  we  see  that  we  must  invert  I  —  AtMA1^ - At2M-1^-.  The  second  term,  —  AfMA1^  = 

— AtM-1D,  is  symmetric  positive  semidefinite  under  a  mass  weighted  inner  product,  which  can  be  used  for 
solving  the  equation.  The  third  term,  however,  can  cause  problems.  Note  that  ^  The  first 

term  is  symmetric  since  Fe  =  —  ,  but  it  might  be  indefinite.  The  second  term  need  not  even  be  symmetric 

under  the  mass  weighted  inner  product. 

Another  approach  is  a  Newmark-style  semi-implicit  scheme 

•  v?+1/2  =  v11  +  ^M-1F(x",vr+1/2) 

•  x"+1  =  xn  +  A  tv?+1/2 

•  v"+1/2  =  vn  +  ^M~1F(x”,vn) 

•  v"+1  =  vn+1/2  +  ^M_1F(x"+1,  v”+1) 
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Each  time  step  of  this  scheme  involves  two  implicit  solves,  that  is  bulleted  steps  1  and  4.  These  solves  are 
different  from  the  fully  implicit  backward  Euler  solve  in  that  the  force  dependence  on  position  is  explicit. 
The  first  solve,  if  taken  for  a  full  At,  would  look  like 


n+1 


n+ 1 


I  -  AtM 


_-yd¥r 


n+1 


<9v 

(I  —  AtM_1D)v”+1 


vn  + AtM-1F(xn,vn+1) 

pn 

vn  +  AtM-1(Fn  +  — —  (vn+1  -  v")) 

av 

^-pn 

vn  -  AtM-1— — v”  +  AtM-1F” 
av 

vn  -  AtM-1Dv"  +  AtM-1F" 


(6) 

(7) 


We  see  from  this  that  the  fully  implicit  system  Eq.  (5)  differs  from  the  semi-implicit  system  Eq.  (6)  by  the 
presence  of  the  extra  —  At2M-1^-  term.  This  term  was  the  source  of  the  indefiniteness  and  asymmetry 
problems  in  the  fully  coupled  system.  This  simplification  comes  at  the  price  of  stability,  however.  The  fully 
implicit  backward  Euler  scheme  is  unconditionally  stable,  but  the  semi-implicit  scheme  is  not.  Of  course, 
this  is  the  well-known  difference  between  a  fully  implicit  Newmark  method  and  one  that  resembles  central 
differencing  [18]. 

The  particular  scheme  we  used  for  the  simulations  in  this  paper  is  based  on  the  backward  Euler  variant 
of  [28],  which  differs  from  the  previous  scheme  in  that  it  has  been  rewritten  so  that  a  backward  Euler  step 
can  be  used  for  updating  the  final  positions.  That  is 

1.  v*  =  v"  +  ^M-1Fe(xn) 

2.  vn+1/2  =  v*  +  +M-1D(xn)v"+1/2 

3.  x”+1  =  x”  +  Afv"+1/2 

4.  v*  =  v"  +  AtM-1Fe(xn+1) 

5.  vn+1  =  v*  +  AtM-1D(x”+1)v”+1 

As  before,  both  solves  lead  to  systems  of  the  form  Eq.  (7).  This  scheme  is  first  order  accurate  in  time,  limited 
to  first  order  by  the  backward  Euler  used  to  compute  the  final  velocities. 

Our  treatment  of  rigid  bodies  follows  that  of  [28].  In  general,  each  particle  of  a  deformable  body  con¬ 
tributes  one  block  to  v,  M,  and  F.  Each  rigid  body  contributes  two  blocks  to  each,  with  the  first  block 
corresponding  to  the  linear  part,  v, ,  m,,  and  fj,  and  the  second  block  corresponding  to  the  angular  part,  ut, 
I i,  and  T*.  In  particular,  the  linear  system  Eq.  (7)  is  valid  for  deformable  and  rigid  bodies. 

The  ability  to  treat  rigid  bodies  using  the  same  formulas  as  deformable  bodies  breaks  down  somewhat 
when  dealing  with  things  like  advancing  rigid  body  orientation.  The  coupling  formulation  presented  here  is 
explicit  in  position  and  therefore  does  not  directly  involve  the  integration  of  rigid  body  orientation,  so  this 
complication  is  avoided.  The  coupling  formulation  is  valid  when  solid  objects  are  deformable,  rigid,  or  a 
mixture  of  both. 


4.  Discrete  Fluid  Structure  Coupling  Equations 

Our  goal  is  to  write  down  the  equations  that  need  to  be  solved  in  full,  including  the  coupling  and  any  other 
terms  such  as  boundary  conditions.  In  doing  this  it  is  more  useful  to  work  from  the  actual  discretizations 
than  it  is  from  the  continuous  equations.  Therefore  we  will  consider  the  MAC  grid  fluid  discretizations  and 
the  Lagrangian  structure  discretizations  as  defined  in  Sections  2  and  3  and  use  them  going  forward  for  our 
notation,  but  we  stress  that  one  could  use  other  discretizations  for  the  fluids  and  structures  and  apply  the 
same  methodology  that  we  propose  to  those  sets  of  equations. 

We  follow  a  method  similar  to  that  of  [25]  to  merge  the  time  integration  schemes  described  in  Sections  2 
and  3, 

1.  v*  =  v"  +  ^M-1Fe(xn) 

2.  u*  =  u”  +  +f 

zp 
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3.  Coupled  solve  with  known  quantities  u*,  v*,  and  x"  and  unknown  quantities  un+1/2,  p,  and  vn+1/2 

4.  xn+1  =  x"  +  Atv"+1/2 

5.  v*  =  v"  +  AtM_1Fe(xn+1) 

6.  u*  =  un  -  At(iin  •  V)u”  + 

7.  Coupled  solve  with  known  quantities  u*,  v*,  and  xn+1  and  unknown  quantities  un+1,  p,  and  vn+1 

Steps  1  through  2  integrate  the  explicit  structure  forces  and  explicit  fluid  forces  for  At/2.  We  leave  out 
the  convective  terms  in  the  fluid  update  because  the  coupled  solve  in  step  3  takes  positions  at  time  n  and 
solves  for  velocities  at  time  n  +  4-  The  known  quantities  at  the  beginning  of  step  3,  which  replaces  step  2 
in  Section  3,  are  u*,  v*,  and  x”.  The  remaining  unknown  quantities  to  solve  for  in  Eqs.  (2),  (4)  and  (7)  are 
Un+i/2,  p,  and  v"+1/2.  We  will  discuss  how  the  actual  coupled  solve  is  performed  later,  except  that  this 
solve  is  performed  over  a  step  size  of  In  step  4,  vra+1/2  is  used  to  advance  the  structure  positions  to  time 
n+  1,  and  then  u™+1/2,  p  and  v”+1/2  are  discarded. 

Steps  5  through  6  re-integrate  the  explicit  structure  forces  and  explicit  fluid  forces  for  At,  this  time 
including  the  convective  terms  for  the  fluid  update  to  get  time  n  +  1  fluid  positions,  resulting  in  v*  and  u*. 
The  known  quantities  at  the  beginning  of  step  7,  which  replaces  step  5  from  Section  3  and  Eq.  (4),  are  u*, 
v*,  and  xn+1.  We  use  x”+1  to  determine  the  valid  degrees  of  freedom  for  u”+1.  The  remaining  unknown 
quantities  to  solve  for  in  Eqs.  (2),  (4)  and  (7)  are  u”+1,  p,  and  vn+1.  The  coupled  solve  used  here  has  the 
same  form  as  that  in  step  3. 

For  the  sake  of  exposition,  we  consider  below  only  the  second  of  the  two  implicit  steps  described  above 
(step  7),  and  therefore  we  are  going  to  discuss  the  updates  from  u*  and  v*  to  un+1  and  vn+1,  whereas  the 
updates  from  u*  and  v*  to  u,i+1/2  and  vn+1/2  will  be  handled  independently.  One  can  replace  At,  x”+1, 
u*,  v*,  u”+1  and  vn+1  with  x”,  u*,  v*,  un+1/2  and  vn+1/2  in  order  to  apply  the  same  treatment  to 
step  3. 

4- 1 ■  Fluid  Structure  Coupling  Constraints 

The  modifications  to  the  fluid  and  structure  equations  correspond  to  constraints  that  capture  the  behavior 
at  the  fluid  structure  interface.  A  fluid  structure  coupling  constraint  is  an  enforced  constraint  on  the  relative 
velocity  of  the  fluid  and  the  structure  at  a  point  in  space.  In  our  implementation,  a  constraint  occurs  when  a 
ray  cast  between  two  adjacent  MAC  grid  cell  centers,  at  least  one  of  which  is  in  the  fluid  region,  intersects  a 
structure  surface.  The  constraint  is  placed  at  the  center  of  the  MAC  grid  face  between  the  two  cell  centers, 
and  enforces  that  the  projection  of  the  fluid  and  structure  velocities  at  that  point  in  the  ray  direction  are 
equal.  The  impulse  transfer  between  the  solid  and  the  fluid  which  enforces  the  constraint  is  called  A.  See 
Figure  1.  This  simple  choice  of  constraint  location  is  only  first  order  accurate,  since  the  constraint  locations 
are  located  0{dx)  from  the  structure’s  surface. 

4-2.  Modified  Fluid  Equations 

For  fluid  structure  interaction,  Eq.  (4)  changes  due  to  the  interaction  with  the  structure.  In  Figure  1  for 
the  x  direction  we  have  a  mass  in  that  dual  cell  in  and  a  velocity  u,  the  product  of  which  mu  is  equal  to 
mu*  plus  the  contribution  from  p  in  the  fluid  on  the  left  side  of  the  dual  cell,  plus  a  contribution  from  A  due 
to  interaction  with  the  structure.  This  discards  the  pressure  on  the  right  hand  side  of  the  dual  cell  due  to 
the  structure.  We  can  neglect  the  dual  cells  that  would  lose  two  pressures,  those  inside  the  structure,  since 
they  are  not  part  of  our  degrees  of  freedom.  We  write  this  in  vector  form  for  all  of  the  equations  as 

(3un+1  =  f3u*  -  Gp  +  WtA,  (8) 

where  (3  is  the  diagonal  matrix  of  fluid  dual  cell  masses  assembled  from  the  product  of  the  density  p  and  the 
volume  V  of  each  dual  cell.  We  let  p  =  Atp.  We  define  — GT  =  —WT,  the  volume  weighted  divergence, 
with  the  additional  modification  that  we  drop  rows  from  —  GT  that  correspond  to  cells  without  fluid.  Its 
negative  transpose,  G,  is  defined  as  the  volume  weighted  gradient.  It  is  simpler  to  describe  —  GT  directly 
rather  than  G,  since  every  row  of  —  GT  has  the  same  stencil  while  the  stencils  for  rows  of  G  differ  due  to 
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Figure  1:  A  dual  cell  requiring  an  equal  velocity  constraint  is  depicted  above  intersecting  the  solid  depicted  by  the  dashed  line. 
The  fluid  pressure  on  the  left  face  of  the  dual  cell  will  be  balanced  by  the  constraint-maintaining  impulse  A  applied  at  the  equal 
velocity  constraint  location. 


some  pressure  samples  falling  inside  the  structure.  WT  is  the  matrix  of  Is  and  Os  which  maps  A  to  the 
appropriate  fluid  velocity  scalars,  which  will  be  discussed  in  more  detail  below.  The  discrete  form  of  the 
incompressibility  constraint  is 

—  Gtu"+1  =  0  (9) 

Our  discretizations  of  W,  G,  and  — GT  are  first  order  accurate  since  the  constraint  locations  used  in 
their  definition  are  only  first  order  accurate.  We  note  that  G  being  a  higher  order  accurate  approximation  of 
the  gradient  operator  does  not  necessarily  imply  that  —  GT  will  be  a  compatible  and  higher  order  accurate 
approximation  of  the  divergence.  Unlike  the  second  order  discretizations  [11,  22],  it  is  not  sufficient  in  this 
context  to  discretize  the  Poisson  operator  with  higher  order  accuracy.  The  discretizations  of  gradient  in 
divergence  in  [22]  do  not  satisfy  the  negative  transpose  relationship,  so  using  them  in  this  formulation  would 
produce  an  asymmetric  linear  system. 

4-3.  Modified  Structure  Equations 

To  write  the  modified  equations  for  the  structure,  one  needs  to  consider  the  fact  that  the  fluid  and 
structure  degrees  of  freedom  are  defined  at  different  spatial  locations.  This  is  quite  important  as  constraints 
must  be  applied  to  co-located  fluid  and  structure  velocities,  and  the  equal  and  opposite  forces  applied  to  the 
fluid  and  the  structure  must  be  applied  at  the  same  location  in  space  so  as  not  to  introduce  a  net  torque. 
We  define  an  operator  J  which  interpolates  from  structure  degrees  of  freedom  to  fluid  velocity  locations.  Its 
transpose  JT  is  a  conservative  redistribution  from  fluid  velocity  locations  to  structure  velocity  degrees  of 
freedom. 

J  can  be  decomposed  into  three  matrices,  J  =  NJR.  The  matrix  R  is  comprised  of  blocks  Rph,  where 
body  b  is  sampled  by  particle  p.  For  particles  within  a  deforming  body,  the  particles  are  the  samples,  and 
Rpb  =  1.  If  b  corresponds  to  a  rigid  body,  then  there  is  a  sample  point  p  at  each  nearby  fluid  velocity 
location.  If  such  a  sample  point  p  is  displaced  from  the  rigid  body’s  center  of  mass  by  rpb,  the  velocity  at 
the  sample  point  is 


vp  =  v6  +r*lub 


(I 


b 


so  that  Rpb  =  (I  r *pb)1  where  r*  represents  the  matrix  such  that  r*u  =  r  x  u  for  any  u. 


(10) 
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The  operator  J  interpolates  particle  velocity  samples  p  to  a  MAC  grid  face  /  as  5^Pa/pvp>  so  that 
the  blocks  of  J  are  J/p  =  a/p I.  Since  J  is  an  interpolation  operator,  J2Pafp  =  1-  The  coefficients  a/p  are 
defined  based  on  local  geometric  information  associated  with  the  face  /  and  the  particle  p  such  as  the  control 
volume  around  /  and  structure  faces  incident  on  p.  For  the  rigid  body,  R  is  used  to  sample  the  rigid  body 
velocities  at  sample  points,  and  the  most  logical  positions  for  sample  points  are  the  fluid  velocity  locations 
at  the  centers  of  faces.  In  that  sense,  if  only  one  rigid  body  was  present  near  a  face,  locally  the  J  operator 
is  the  identity  for  that  degree  of  freedom.  However  multiple  rigid  bodies  or  particles  may  influence  the  same 
face,  in  which  case  an  averaging  of  the  velocities  would  be  required,  making  J  nontrivial. 

Finally,  N  extracts  the  appropriate  scalar  velocity  component  from  the  vector  velocity  constructed  via 
R  and  J.  The  matrix  N  is  block  diagonal,  where  each  block  N ff  corresponds  to  a  grid  face  /  with  normal 
ii/.  These  blocks  are  formed  as  N //  =  nj.  Our  discretization  of  the  operator  J  is  limited  to  first  order 
accuracy  since  it  depends  on  the  first  order  accurate  structure  interface  location. 

Since  Eq.  (8)  applies  an  impulse  to  the  fluid  at  MAC  grid  face  locations,  the  equal  and  opposite  impulse 
must  be  applied  to  the  structure  at  those  same  locations,  guaranteeing  conservation  of  angular  momentum. 
For  the  case  of  rigid  bodies,  one  can  apply  impulses  to  the  rigid  body  from  any  location  in  space  and  therefore 
the  method  does  conserve  angular  momentum,  however  in  the  case  of  deformable  bodies  JT  contains  a  spatial 
redistribution  of  the  impulse  which  while  being  a  fully  conservative  remapping  with  columns  summing  to  1 
and  thus  conserving  linear  momentum,  will  not  conserve  angular  momentum.  In  summary,  NT  maps  impulse 
magnitudes  from  dual  cells  to  impulse  vectors  normal  to  the  face.  Since  J2pctfp  =  1,  JT  conservatively 
distributes  these  impulses  to  different  spatial  sample  locations  in  the  case  of  deformable  bodies.  For  rigid 
bodies  the  spatial  location  does  not  change  although  the  impulse  can  still  be  conservatively  distributed 
among  multiple  rigid  bodies  influencing  that  face.  Finally  RT  conservatively  applies  those  sample  impulses 
to  the  structure  degrees  of  freedom.  This  mapping  is  trivial  for  deformable  particles,  but  for  rigid  bodies  it 
converts  the  impulse  sample  into  linear  and  angular  impulses  at  the  body’s  center  of  mass.  We  note  that  [2] 
proposed  a  method  which  allows  the  solid  and  fluid  domains  to  be  non-coincident,  coupled  via  interpolation 
and  Lagrange  multipliers,  similar  t  our  use  of  W  and  J. 

The  full  modified  structure  equations  are 

Mv"+1  =  Mv*  +  AtDv"+1  —  JtWtA  (11) 

where  M  is  the  structure  mass  matrix  and  Dvn+1  is  the  implicit  contribution  of  structure  forces.  The 
last  term  is  the  constraint  impulse  mapped  first  by  WT  from  constraint  locations  to  fluid  velocity  degrees 
of  freedom  and  then  by  JT  to  structure  velocity  degrees  of  freedom.  The  coupling  momentum  WrA  is 
subtracted  from  the  structure  equation  since  it  was  added  to  the  fluids  equation.  W  has  a  number  of 
columns  equal  to  the  number  of  fluid  scalar  velocity  samples  and  a  number  of  rows  equal  to  the  number  of 
constraints.  WT  therefore  maps  to  every  MAC  grid  face  but  puts  zero  on  every  face  which  is  not  coupled. 
As  a  consequence,  JT  only  needs  to  be  nonzero  for  the  faces  which  are  coupled,  even  though  it  has  a  number 
of  columns  equal  to  the  number  of  fluid  scalar  velocity  samples.  JTWT  then  maps  from  the  fluid  degrees  of 
freedom  that  are  coupled  to  the  structure  degrees  of  freedom  that  are  coupled. 

Although  Eqs.  (8)  and  (11)  indicate  how  momentum  is  conservatively  transferred  between  the  fluid  and 
the  structure,  this  still  leaves  open  the  question  of  how  much  momentum  is  transferred.  That  is,  we  have 
introduced  a  new  degree  of  freedom  A  and  have  not  yet  provided  the  equations  to  govern  these  new  degrees 
of  freedom.  In  order  to  do  this  we  consider  the  equal  velocity  constraints  that  the  fluid  and  the  structure 
should  move  with  the  same  velocity,  or  for  inviscid  flow  the  same  velocity  in  the  normal  direction.  One  way 
of  governing  these  degrees  of  freedom  is  to  set  the  difference  between  the  fluid  velocity  and  the  interpolated 
structure  velocity  equal  to  zero  at  every  fluid  velocity  degree  of  freedom  as  indicated  in  Figure  1.  This  can 
be  done  in  equation  form  via 


W(u-Jv)  =  0,  (12) 

noting  that  J  maps  zero  to  all  the  face  locations  that  are  not  influenced  by  the  structure  and  that  W  filters 
all  of  these  locations  out.  Note  that  one  can  use  Eq.  (12)  and  its  variants  in  order  to  model  a  wide  range 
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of  physical  phenomena.  For  example,  one  could  use  W(u  —  J(v  —  s))  =  0  in  order  to  model  an  arbitrary 
source  term  with  inflow  and  outflow  from  the  structure.  One  might  also  wish  to  constrain  the  velocities 
at  other  locations  than  the  faces  of  the  MAC  grid  or  to  constraint  limited  components  such  as  the  normal 
component  as  was  done  in  [24]  in  order  to  achieve  better  slip  boundary  conditions.  In  this  sense  W  may  not 
be  simply  defined  as  a  trivial  restriction  matrix  and  may  instead  be  thought  of  as  a  general  interpolation 
operator  mapping  from  the  faces  of  the  MAC  grid  to  constraint  locations.  This  would  not  change  the  inputs 
of  W,  but  the  outputs  would  no  longer  necessarily  correspond  to  the  MAC  grid  faces  but  could  instead  be 
at  arbitrary  locations.  Note  that  the  formulation  of  our  equations  does  not  restrict  any  of  these  approaches, 
but  for  the  sake  of  exposition  we  will  deal  with  Eq.  (12)  going  forward. 

5.  Numerical  Approach 

We  can  write  our  system  in  symmetric  form  as 

/  GT/3-1G  -Gt^_1Wt  0  \  /  p  \  /  Gtu*  \ 

— W/3_1G  W/3_1Wt  -WJ  A  =  -Wu*  . 

V  0  -JTWT  — M  +  At  D/  \vn+7  \-M V*/ 

The  first  equation, 

GT/3~1Gp-GT/3"1WTA  =  GTu*  (13) 

is  obtained  by  substituting  Eq.  (8)  into  Eq.  (9).  The  second  equation, 

— W/3_1Gp  +  W/3~1WtA  -  WJv"+1  =  -  Wu*  (14) 

is  obtained  by  substituting  Eq.  (8)  into  Eq.  (12).  The  third  equation, 

-JtWtA  -  Mv"+1  +  Af  Dv"+1  =  -Mv*  (15) 

is  a  rearrangement  of  Eq.  (11). 

If  we  make  some  assumptions  about  W  and  J,  we  can  eliminate  A.  If  we  assume  the  W  takes  on  the 
particularly  simple  form  that  we  used,  then  WWT  =  I.  If  J  has  nonzero  rows  only  for  constraint  faces,  then 
WTWJ  =  J.  If  we  also  take  densities  to  be  constant  at  all  constraint  faces,  then  WTW/3_1  =  /3_1WTW. 
Multiplying  Eq.  (14)  by  WT  produces 

-WTW/3"1Gp  +  WTW/3~1WTA- WTWJvn+1  =  WTWu* 

-WTW/3"1Gp  +  /3"1WTA- Jv"+1  =  WTWu*  (16) 

Adding  GT  times  Eq.  (16)  to  Eq.  (13)  yields 

(Gt/3-1G  —  GTWTW/3-1G)p  +  (— Gt/3_1Wt  +  GT/3-1Wr)A  —  GTJv”+1  =  Gru*  -  GTWTWu* 

Gt(I  -  WTW)/3_1Gp  -  GtJv”+1  =  Gt(I  —  WTW)u* 

Adding  JT/3  times  Eq.  (16)  to  Eq.  (15)  yields 

— JT/3WTW/3_1Gp  +  (-JTWT  +  JTWT) A  +  (At  D  -  M  -  JT/3J)v"+1  =  -Mv*  +  JT/3WTWu* 

-JTGp  +  (AfD-M-  Jt/3J)v"+1  =  -Mv*+Jt/3u* 

These  together  form  a  symmetric  system  similar  to  that  of  [25] 

/Gt(I-  WtW)/3"1G  -GtJ  )  /  p  \  /GT(I- WTW)u*\ 

V  JTGp  -Jt/3J-M  +  AtD/  \vn+1)  ~  V  -Mv*  +  JT/3u*  J  ' 
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Both  of  these  system  are  symmetric  indefinite.  Although  numerical  methods  exist  to  solve  these  systems, 
faster  methods  are  available  for  symmetric  positive  definite  systems.  In  this  paper,  we  propose  a  novel 
approach  which  leads  to  a  symmetric  positive  definite  system. 

We  begin  by  forming  the  system 

/  Gt(3~1G  —Gt (3~1Wt  0  \  /  p  \  /  Gtu*  \ 

-W/3_1G  W^V+WJM'^W1,  — AtWJM_1D  I  I  A  =  WJv*  -  Wu*  (17) 
V  0  M^1JtWt  I  —  At  M-1D  /  \vn+1J  V  v*  / 

where  the  second  equation,  -W/3_1Gp  +  W(3~1WT\  +  WJM^J^^- At  WJM-'Dv^1  =  WJv*- 
Wu*,  is  obtained  by  substituting  Eq.  (8)  and  Eq.  (11)  into  Eq.  (12).  The  first  equation  is  Eq.  (13),  and  the 
third  equation  is  Eq.  (15)  multiplied  through  by  M-1. 

An  obvious  way  to  symmetrize  this  system  is  to  multiply  the  third  equation  through  by  AID.  Though 
this  system  is  symmetric  it  is  only  positive  semidefinite.  The  conditioning  of  the  original  structure  system  is 
typically  quite  good,  since  I  —  At  M'1D  has  eigenvalues  clustered  near  one  for  small  At.  This  symmetrized 
version  looks  more  like  D,  which  will  often  have  significantly  worse  conditioning. 

A  particularly  interesting  way  of  symmetrizing  the  system  is  to  factor  D. 

5.1.  Factoring  Damping  Matrix 

If  D  is  symmetric  negative  semidefinite,  then  it  can  be  factored  as  —  AtD  =  CTC.  We  rely  on  this 
property  of  D  to  apply  our  method,  noting  that  although  there  are  instances  where  D  will  not  have  this 
property  and  in  those  cases  our  approach  based  on  factoring  D  is  not  applicable,  there  are  a  wide  range 
of  models  in  which  this  property  is  true  (see  below).  If  D  is  a  damping  matrix,  and  the  force  is  given  by 
f  =  Ma  =  —  \7<f>  +  Dv,  then  the  energy  is  E  =  tp  +  fvTMv  and  E  =  (V</>)Tv  +  vTMa  =  ( V(f>)Tv  + 
vT(-V0  +  Dv)  =  vtDv.  Then,  a  non-positive  change  in  energy  leads  to  negative  definiteness  of  D. 
Furthermore,  if  D  is  decomposed  into  its  symmetric  and  antisymmetric  parts,  D+2D  and  D  ,,D  respectively, 
then  vt^d~2d  ^jv  =  0  shows  that  the  antisymmetric  part  does  not  contribute  damping.  Therefore,  damping 
matrices  gain  little  by  containing  asymmetric  parts.  These  facts  lead  to  the  conclusion  that  for  most  purposes 
a  damping  matrix  D  can  be  considered  to  be  symmetric  negative  semidefinite. 

On  the  other  hand,  if  you  consider  the  backward  Euler  formulation  at  the  beginning  of  Section  3,  where 
D  would  be  augmented  by  linearizations  of  elastic  terms  in  the  process  of  carrying  out  the  Newton- 
Raphson  iterations,  then  much  of  this  would  not  be  true.  The  current  factorization  is  therefore  limited  to  a 
semi-implicit  treatment  of  the  structure,  where  elastic  terms  are  handled  explicitly.  Moreover,  a  fully  implicit 
treatment  would  have  the  further  complication  of  moving  the  locations  of  the  constraints  and  changing  W 
during  the  solve,  which  we  do  not  address. 

For  the  sake  of  exposition,  consider  the  case  of  a  linear  spring  which  connects  two  particles  pi  and  P2. 
The  magnitude  of  the  damping  force  f  is  proportional  to  the  difference  in  particle  velocities  along  the  spring 
direction  d,  multiplied  by  a  damping  constant  kd ■  The  force  is  then  applied  in  opposite  directions  to  pi  and 
P2 .  The  damping  matrix  for  this  single  spring  would  be 


which  can  be  factored  into 


Cj  = 

D  =  -CfCi. 

For  the  case  of  n  springs  and  an  arbitrary  number  of  particles,  C  is  formed  by  stacking  up  rows  of  the  form 
C i  and  the  same  factorization  applies. 
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Rayleigh  damping  is  commonly  used  along  with  finite  element  constitutive  models.  In  practice,  we  have 
found  that  our  routines  for  applying  damping  forces  are  naturally  partially  factored.  The  implementation 
follows  the  framework  of  [20].  Below  is  a  description  of  what  the  implementation  looks  like  in  the  case  of  a 
single  element  and  an  isotropic  constitutive  model  and  how  it  lends  itself  to  efficient  factorization.  In  this 
framework,  the  deformation  gradient  is  factored  into  its  polar  singular  value  decomposition  as  F  =  UFVT, 
where  U  and  V  are  rotations  and  F  is  diagonal,  in  order  to  compute  the  elastic  forces.  Define  two  operators 


=  (xi  -  x4  x2  —  x4  x3  -  x4)  S  (xi  x2  x3)  = 

If  xm  are  the  material  space  positions  of  the  element,  then  Dm  =  G(xrn)  and  Ds  =  G(x)  are  the  locations 
of  three  vertices  of  the  element  relative  to  the  fourth  in  material  and  world  space.  The  damping  force 
application  routine  consists  of  these  steps 

1.  Compute  the  relative  velocities,  Ds  =  G(v). 

2.  Compute  the  deformation  gradient,  F  =  DgD”1. 

3.  Rotate  the  deformation  gradient  into  the  diagonalized  coordinate  system,  F  =  UTFV. 

4.  Compute  the  first  Piola-Kirchhoff  stress  in  the  diagonalized  coordinate  system,  P  =  atr(F)I  +  (3F  + 

(3FT ,  where  a  and  (3  are  parameters  of  the  constitutive  model. 

5.  Rotate  back  to  compute  the  first  Piola-Kirchhoff  stress,  P  =  UPVT. 

6.  Compute  forces  the  element  applies  on  three  of  the  tetrahedron  faces  as  G  =  —  PN,  where  N  are 
the  area  weighted  face  normals.  These  normals  can  be  computed  as  N  =  —  ^  det(Dm)D“T.  To  see 
that  this  is  correct,  note  that  |NTDm  =  —VI,  where  V  is  the  volume  of  the  tetrahedron.  The 
negative  sign  is  because  N  need  to  point  outwards.  Note  that,  if  we  computed  it,  G4  =  -PN4  = 
-P(-N1-N2-N3)  =  -Gi-G2-G3. 

7.  Distribute  the  forces  from  faces  to  vertices.  The  force  f4  at  the  first  vertex  is  one  third  of  the  forces  on 
the  faces  around  it,  or  f4  =  i(G2  +  G3  +  G4)  =  —  |Gi-  Noting  that  the  forces  sum  to  one,  the  forces 
are  finally  obtained  as  f  =  S(—  |G). 

Several  observations  make  the  process  of  factoring  this  force  application  simpler.  The  first  observation  is 
that  each  of  these  steps  can  be  regarded  as  linear  operators  between  vector  spaces,  where  the  matrices  are 
taken  to  be  flattened  into  vectors.  In  this  way,  the  operators  S  and  G  are  transposes,  so  that  ignoring  the 
factor  of  —  the  operator  described  by  step  7  is  the  transpose  of  step  1.  Similarly,  step  6  is  \  det(Dm)  times 
the  transpose  of  step  2.  Step  3  is  the  transpose  of  step  5.  Let  P{a,  (3)  represent  the  linear  operator  applied  in 
step  4.  The  damping  force  application  can  be  written  as  steps  1  through  3,  followed  by  —  \  det(Dm)P(a,  (3), 
followed  by  the  transpose  of  steps  1  through  3.  Finally,  the  operation  C  is  just  steps  1  through  3,  followed 
by  the  square  root  of  ^  At  det(Dm)P(a,  (3).  It  remains  only  to  determine  what  this  square  root  is. 

Assume  that  the  square  root  has  the  form  P(a'  ,(3').  Then, 

P(a',f3')P  =  a'tr(^)I  +  +  (3'tT 

P(a  ,  (3')2 F  =  c/Itr(c/tr(F)I  +  /3' F  +  (3'Ft)  +  /3'(c/tr(F)I  +  0F  +  /3'Ft)  +  /3'(c/tr(F)I  +  /3'F  +  /3"Ft)t 
=  (3a'2  +  4c//3')tr(F)I  +  2/3'2F  +  2/3'2FT 
=  P(3c/2  +  4a'f3',  2(3,2)F 

(3 '  =  ^Atdet(Dro)(iy/2^) 

Q(\/3^+2^- 

Since  the  operator  P(a',(3')  is  symmetric,  this  construction  results  in  the  operator  C,  and  routines  to  apply 
both  C  and  CT  are  readily  constructed  with  only  minor  modifications  to  the  existing  code. 
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5.2.  SPD  Formulation 

Assuming  a  factorization  of  the  form  —  At  D  =  CTC  allows  us  to  left  multiply  the  third  row  of  Eq.  (17) 
by  C  and  subsequently  factor  C  out  of  the  third  column  into  the  solution  vector  to  obtain 


/  Gt/3_1G  — Gt/3-1Wt  0 

-W/3_1G  W/3_1Wt  +  WJM-1JtWt  —At  WJM-1D 
V  0  M_1JtWt  I  AtM^'D 


/  GT/3~1G  — GT/3_1WT  0 

-W^G  Wf3~1WT  +  WJM-jJtWt  WJM1CrC 
V  0  M_1JtWt  I  +  M~1CtC 


/  Gt/3~1G 

— Gt/3-1Wt 

°  \  { 

-Wfi^G 

W(/3_1  +  JM-1Jt)Wt 

WJM^C7 

V  0 

cm-1jtwt 

I  +  CM^Ct/  \ 

(  GTu *  \ 
I  WJv*  Wu*  j 


/  GTu*  \ 

I  WJv*  Wu*  j 

(  Gtu*  \ 

WJv*  Wu* 

V  Cv*  j 


where  v  =  Cvn+1  represents  the  new  unknowns.  Note  that  the  substitution  of  Eq.  (11)  into  Eq.  (12)  to 
derive  Eq.  (17)  allowed  us  to  place  a  D  into  the  third  column  of  the  coefficient  matrix  in  Eq.  (17),  leading 
to  the  factorization  of  C  out  of  the  coefficient  matrix  and  into  the  unknowns,  and  this  was  the  reason  for 
that  extra  substitution.  The  resulting  system  is  obviously  symmetric,  and  we  can  recover  v"+1  from  v  by 
replacing  AfDv"+1  with  — CTCv"+1  =  — CTv  in  Eq.  (11) 

Mv"+1  =  Mv*  -  CTv  -  JtWtA. 


The  symmetrized  system  can  be  expressed  as  the  sum  of  three  matrices  that  are  evidently  symmetric 
positive  semidefinite  as 

/GTf3~1G  -Gt/3~1Wt  0\  /0  0  0  \  /0  0  0\ 

— W/3-1G  W/3-1Wt  0  +  0  WJM  !JTWT  WJM^C1"  +  0  0  0 

\  0  0  0/  V°  CM-^W7  CM-'Ct  J  \0  0  I ) 

from  which  the  symmetric  positive  semidefiniteness  of  the  coupled  system  follows. 

To  consider  the  rank  of  this  matrix,  note  that  any  vector  in  the  nullspace  of  the  coupled  system,  consisting 
of  three  parts  p,  A,  and  v  must  necessarily  be  in  the  nullspace  of  each  of  these  three  matrices.  From  the 
third  matrix,  we  conclude  that  this  requires  v  =  0.  Assuming  v  =  0  and  looking  at  the  second  matrix  leads 
one  to  only  consider  requirements  for  A.  The  second  matrix  factors  as 

°  \  /  0  \T 

WJ  M  1  WJ 

c  )  \c  J 

It  is  necessary  for  M-1JTWTA  =  0,  which  implies  that  the  constraint  impulses  have  absolutely  no  effect  on 
any  degrees  of  freedom  of  the  structure.  The  first  matrix  factors  as 

GT\  f Gt\T 

W  /3-1  W 

0  /  V  0  J 

A  nontrivial  nullspace  of  this  matrix  satisfies  /3-1Gp  /3_1WTA  =  0,  which  implies  that  the  combination  of 

constraints  and  pressures  cancels  to  give  no  effect  on  the  final  fluid  velocities  u"+1  as  can  be  seen  by  inspecting 
Eq.  (8).  This  immediately  gives  the  cases  /3_1Gp  =  0  (Neumann  region)  and  WTA  =  0  (redundant 
constraints),  both  of  which  will  lead  to  nullspaces  of  the  full  coupled  system  by  placing  zeros  in  the  remaining 
entries. 
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The  right  hand  side  of  our  equations  must  be  in  the  range  of  the  coefficient  matrix  in  order  for  the 
system  to  be  solvable.  Comparing  the  form  of  our  right  hand  side  to  the  requirements  that  a  nullspace  must 
satisfy,  we  find  that  it  is  sufficient  but  not  necessary  condition  for  consistency  for  v*  to  contain  no  nullspace 
components  of  MA1  and  u*  to  contain  no  nullspace  components  of  /3-1.  A  Neumann  region  will  result  in 
inconsistency  if  there  is  a  net  flux  in  u*  through  the  region  boundary,  so  we  enforce  consistency  by  summing 
up  the  net  flux  and  subtracting  it  evenly  from  the  u*  on  the  region  boundary  faces.  Clearly  constraints 
which  apply  only  to  nullspace  components  of  M_1  and  /3_1  will  result  in  inconsistency  and  must  be  avoided. 
Fluid  velocities  which  are  constrained  only  to  nullspace  components  of  MA1  are  in  effect  Neumann  grid  faces 
and  must  be  taken  into  account  in  detection  of  Neumann  regions.  Isolated  rigid  bodies  within  a  fluid  region 
also  do  not  prevent  a  region  from  being  a  Neumann  region.  A  region  that  contains  deformable  bodies  does 
not  behave  as  a  Neumann  region,  since  the  bodies  can  change  volume  to  accommodate  a  net  flux  into  the 
region. 


5.3.  Note  On  Forming  SPD  System 

Many  problems  can  be  formulated  as  a  KKT  system  of  the  form 


(M  +  CTC  Bt\  /u\  _  /b\ 

l  B  0  )  [x)  ~  [o) 


(18) 


where  M  is  SPD  matrix  that  is  easily  inverted.  Dissipative  systems,  such  as  fluid  with  viscosity  or  structure 
damping,  produce  terms  of  the  form  M  +  CTC.  Evolving  such  systems  subject  to  constraints,  such  as  the 
incompressibility  constraint  or  the  constraint  that  the  solid  and  fluid  velocities  be  equal  at  their  interface, 
then  leads  to  a  KKT  system  in  the  form  of  Eq.  (18).  Multiplying  the  first  row  of  Eq.  (18)  by  CM  1  produces 

Cu  +  CM  1CTCu  +  CM  2BtA  =  CM-1b.  (19) 

Multiplying  the  first  row  of  Eq.  (18)  by  BM^1  and  simplifying  using  the  second  row  produces 

BM  1CTCu  +  BM_1BtA  =  BM^b.  (20) 


These  are  then  combined  into  the  SPD  system  using  the  definition  u  =  Cu  to  produce 

/I+CM^Ct  CM^1Bt\  /u\  /CM-1^ 

^  bm-'ct  bm^b7]  va/  vBM_lb/  ' 


Finally,  u  is  obtained  from  the  first  row  of  Eq.  (18)  using 

u  =  M~1b-M  1Ctu-M  1BtA.  (22) 

To  demonstrate  that  this  is  in  fact  equivalent  to  Eq.  (18),  we  show  that  we  can  derive  Eq.  (18)  from  u  =  Cu, 
Eq.  (21)  and  Eq.  (22).  The  first  row  of  Eq.  (18)  is  obtained  by  substituting  u  =  Cu  into  Eq.  (22).  The 
constraint  row  is  obtained  by  subtracting  the  second  row  of  Eq.  (21)  from  B  times  Eq.  (22).  Thus,  the 
indefinite  KKT  formulation  is  equivalent  to  the  SPD  formulation. 


6.  Note  On  Conditioning 

In  the  absence  of  coupling  (W  =  0),  the  structure  system  looks  like 

vn+r  =  v*  +  AtM-1Dvn+1 
vn+1-  AiM_1Dvn+i;  =  V* 

(I  +  M_1CtC)v"+1  =  v* 

The  modified  structure  system  looks  like 

(I  +  CM-1Ct)v  =  Cv* 
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These  two  systems  have  the  same  eigenvalues.  To  see  this 

(I  +  M_1CTC)v  =  av 
C(I  +  M-1CtC)v  =  Cav 
(I  +  CM_1CT)v  =  av 

so  that  if  v  is  an  eigenvector  with  eigenvalue  a  in  the  original  system,  then  v  is  an  eigenvector  of  the  modified 
system  with  the  same  eigenvalue.  Note  also  that,  though  modifying  the  system  in  this  way  may  introduce 
new  eigenvalues  (since  the  system  may  be  larger),  none  of  these  eigenvalues  will  be  visible  to  a  Krylov  solver. 
This  is  because  the  right  hand  side  is  of  the  form  Cv*  and  application  of  the  matrix  results  in  vectors  that 
are  also  of  this  form.  The  modified  system  may  have  fewer  eigenvalues,  however  (caused  by  rank  deficiency 
of  C).  The  important  point  is  that  this  modification  to  the  structure  system  will  not  adversely  affect  its 
convergence  properties  with  a  Krylov  solver. 


7.  Viscosity 

Applying  viscosity  as  a  separate  step  in  the  time  evolution  does  not  yield  a  fully  coupled  time  integration 
scheme,  since  viscosity  forces  cannot  apply  momentum  to  structures.  This  can  be  overcome  by  adding 
viscosity  directly  into  the  coupled  formulation.  With  this  in  mind,  we  redefine  the  intermediate  velocity 

u*  =  u"  -  Af(un  •  V)u”  +  —  f.  (23) 

P 

Assume  that  the  fluid  viscosity  p  is  constant  and  ignore  the  coupling  term.  Then,  the  fluid  momentum 
equation  becomes 

u*  =  u”  -  Af(un  •  V)u"  +  —  f/9un+1  =  pu*  -  A tVp  +  AlpV2un+1.  (24) 

P 

The  viscous  terms  introduce  a  new  gradient  operator,  which  acts  on  face  velocities.  The  gradient  operator 
itself  is  discretized  using  central  differences.  The  boundary  conditions  at  structures  are  Dirichlet  conditions, 
noting  that  these  velocities  are  constrained  implicitly  using  W  and  remain  degrees  of  freedom  rather  than 
being  moved  to  the  right  hand  side  as  would  normally  be  done.  The  discretization  of  W  used  above  is  not 
suitable  for  viscosity,  however,  since  it  does  not  couple  to  structures  in  the  tangential  direction.  This  is 
resolved  by  locating  faces  adjacent  to  those  in  W  and  adding  them  to  W  as  shown  in  Figure  2.  Domain 
boundary  conditions  can  be  discretized  in  the  usual  way.  For  simplicity  of  implementation,  however,  we 
treat  domain  walls  as  infinitely  massive  structures  and  add  terms  to  W,  J,  and  M-1  accordingly.  This 
discretization  is  limited  to  first  order  accuracy  since  it  uses  the  first  order  accurate  constraint  locations. 

Let  i /KtjjVV  be  the  discretized  gradient  operator  that  acts  on  face  velocities,  where  the  scaling 

factors  are  pulled  into  the  resulting  matrix  to  simplify  the  presentation.  Substituting  this  into  Eq.  (24)  and 
adding  back  the  coupling  term  produces  the  momentum  equation 

f3un+1  =  f3u*  Gp  -  G^G^u”-1"1  +  WtA.  (25) 

Define  u  =  G;,un+1 .  Then,  multiplying  Eq.  (25)  through  by  both  GT/3~1  (using  Eq.  (9))  and  G^/3^1 
produces 

GT/3“1Gp  +  Gt/3-1G^u  —  Gt/3-1WtA  =  GTu*  (26) 

GM/3_1Gp  +  (I  +  GAj/3_1GJ)u  —  GM/3_1WTA  =  G^u*  (27) 

Multiplying  Eq.  (25)  through  by  W/3_1  and  substituting  in  Eq.  (12)  followed  by  Eq.  (11),  —  AfD  =  CTC 
and  Cv”+1  =  v  produces 

W/3_1Gp  +  W/3_1G^u  —  W(/3_1  +  JM_1JT)WTA  —  WJM'1CTv  =  Wu*-WJv*  (28) 
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Figure  2:  The  hashed  region  is  the  rasterized  structure,  and  the  thick  faces  at  its  boundary  are  the  original  coupling  faces  in 
W.  For  viscosity,  the  dashed  faces  are  added  to  W  to  allow  tangential  forces  to  be  exchanged  between  fluid  and  structure. 


Combining  these  with  Eq.  (11)  produces 


/  Gt/3_1G 

G+T'GJ 

-Gt/3-1Wt 

°  ^ 

(p) 

(  GTu*  \ 

G^G 

i  +  g+t^ 

-Gm/3-1Wt 

0 

U 

Gmu* 

-Wp^G 

-w/r'G  l 

W(/3-1  +  JM-1Jt)Wt 

wjm-1ct 

A 

WJv*  -  Wu* 

V  0 

0 

cm-'jtwt 

I+CM^CV 

V  Cv*  / 

This  equation  has  a  structure  very  similar  to  what  was  obtained  without  viscosity,  and  it  can  be  shown  to 
be  symmetric  positive  semidefinite  in  exactly  the  same  way. 


Stokes  Flow.  When  the  Reynolds  number  is  very  small,  the  viscosity  term  dominates  the  inertial  terms. 
Because  viscosity  and  pressure  are  fully  coupled  in  a  single  linear  system,  our  scheme  does  not  experience 
difficulties  at  high  viscosity.  In  this  case,  one  may  also  approximate  the  Navier-Stokes  equations  with  those 
of  Stokes  flow  by  omitting  the  inertial  terms.  We  can  simulate  the  unsteady  Stokes  flow  equations  directly 
in  our  framework  by  omitting  advection  terms  for  the  fluid  (Step  6)  and  the  solid  (Step  4).  This  renders  the 
entire  first  half  of  the  update  unnecessary,  and  our  evolution  scheme  reduces  to 

1.  v*  =  v"  +  AtM~1Fe(xn+1) 

2.  u*  =  u"  + 

p 

3.  Coupled  solve  with  known  quantities  u*,  v*,  and  x"+1  and  unknown  quantities  u"+  ,  p,  and  v"+1 

We  demonstrate  our  ability  to  simulate  with  high  viscosity  and  Stokes  flow  in  Section  9.2. 

As  discussed  in  [25],  the  indefinite  system  of  equations  derived  for  fluid-structure  interaction  has  many 
similarities  to  the  Stokes  equations  as  both  are  dissipative  systems  combined  with  constraints  (see  Sec¬ 
tion  5.3).  Indeed,  we  produce  an  SPD  formulation  for  both.  If  no  structures  are  present  then  W  =  0,  C  =  0 
and  J  =  0,  and  Eq.  (29)  simplifies  to 


( GT/3-1G  GyGj\/p\  =  /  GTu*\ 

^G+t'g  I  +  GM/3-1G^y  \uy  V<+W  ' 


(30) 


This  system  provides  a  SPD  treatment  of  Stokes  flow.  If  convection  is  handled  explicitly,  this  formulation 
also  provides  a  SPD  treatment  of  Navier-Stokes. 


Fully  implicit  Navier  Stokes.  We  could  take  a  fully  implicit  approach  to  the  full  nonlinear  Navier-Stokes 
equations  by  using  a  Newton-Raphson  iteration  method,  where  we  linearize  at  each  Newton  step  to  get  a 
system  of  the  form  in  Eq.  (30).  The  equation  to  be  iterated  is 


u 


n+1  _ 


(i+l) 


n+1 


=  u"-At(u~ 


V)^1 


At,  _ 

—  -Vp 

P 


f), 
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where  u^j"1  is  the  time  tn+1  fluid  velocity  after  the  ith  Newton  iteration  and  we  let  u”^1  =  u”.  The  inner 
iterations  are  performed  by  solving  the  symmetric  and  positive  definite  system  Eq.  (30).  For  each  outer 
iteration,  a  new  system  is  set  up  using  the  new  estimate  of  the  final  velocity  to  recompute  the  advection 
term.  A  similar  approach  could  be  explored  for  treatment  of  fully  implicit  coupled  fluid-structure  evolution. 


8.  Note  on  Stability 


We  have  found  the  proposed  monolithic  implicit  coupling  scheme  described  in  this  paper  to  be  quite 
stable.  We  have  never  observed  stability  problems  due  to  this  simple  choice.  In  the  example  from  Figure  9, 
we  consider  the  structure  discretization  to  be  exact  and  refine  only  the  fluid  discretization.  This  results  in  a 
grid  spacing  much  smaller  than  the  lengths  of  edges  in  the  structure  mesh,  which  has  never  been  observed 
to  cause  stability  problems.  We  also  tried  running  examples  with  a  fixed  resolution  fluid  mesh,  refining 
the  structure  so  that  the  grid  spacing  was  over  30  times  the  lengths  of  the  structure  edges,  again  without 
noticing  any  stability  problems.  We  have  tried  making  structures  a  billion  times  less  dense  than  the  fluid 
and  a  billion  times  more  dense  than  the  fluid,  and  neither  modification  lead  to  stability  problems. 

We  also  note  that  solving  the  coupled  linear  system  and  applying  the  resulting  update  the  velocities  of 
the  structure  and  the  fluid  has  the  property  that  it  never  increases  total  kinetic  energy.  Since  positions  are 
not  changed  during  this  process,  potential  energy  is  unchanged.  Let 


K  = 


(Gt 

G„ 

-W 

V  o 


0  \ 

0 

WJ 

C  J 


/o  0  0  0\ 

(p\ 

AS  0\ 

0  10  0 

u 

(o  m)  p  = 

0  0  0  0 

w=W  z= 

A 

^0  0  0  1^ 

w 

Eq.  (29)  reduces  to  (KB  1KT  +  P)z  =  Kw.  The  new  velocities  are  w  —  B  1KTz.  Finally,  the  change  in 
kinetic  energy  is 


2 KEn+1  -  2KE* 


(w  B  1Ktz)tB(w  —  B  1Ktz)  —  wtBw 
— 2wtKtz  +  (Ktz)tB^1(Ktz) 

— 2zt(KB-1Kt  +  P)z  +  ztKB-1Ktz 
-zt(KB_1Kt  +  2P)z  <  0 


This  does  not  imply  that  we  always  decrease  kinetic  energy,  since  elastic  forces  added  to  the  structures  or 
gravity  added  to  fluids  may  both  increase  kinetic  energy.  In  particular,  we  are  still  explicit  on  elastic  forces, 
and  this  continues  to  impose  a  CFL  restriction.  On  the  other  hand,  if  there  is  no  potential  energy  (e.g.,  rigid 
bodies  without  forces  between  them),  and  the  remaining  steps  in  the  time  integration  scheme  make  a  similar 
guarantee  on  kinetic  energy,  the  resulting  scheme  would  be  unconditionally  stable.  We  also  note  that  if 
there  is  neither  damping  nor  viscosity,  then  P  =  0,  and  the  process  of  obtaining  the  final  velocities  from  the 
intermediate  velocities  w  is  described  by  the  mass-synnnetric  projection  matrix  I  B  1Kr(KB  1KT)_1K, 
which  is  exactly  the  property  that  would  be  desired  in  this  case. 


9.  Examples 

The  SPD  formulation  merely  gives  an  analytically  equivalent  expression  of  the  equations.  What  we 
expect  to  illustrate  in  these  examples  is  that  the  method  is  numerically  stable.  We  examine  several  common 
examples  from  the  fluid  structure  interaction  literature  to  demonstrate  the  convergence  properties  of  our 
method. 

We  choose  our  time  step  size  to  be  the  minimum  of  the  time  step  size  required  for  fluids  alone  and  the 
time  step  size  required  for  structures  alone.  The  fluid  time  step  size  is  computed  using  a  CFL  number  of 
0.9  for  all  of  our  examples.  Except  in  the  stability  test  described  below,  the  structures  in  our  examples  are 
not  refined,  so  the  time  step  size  computed  for  the  structures  will  not  change  significantly  from  resolution 
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Figure  3:  Analytic  test  with  a  structure  flowing  through  a  channel.  The  errors  are  computed  as  the  difference  between  the 
steady  state  velocity  of  the  moving  structure  and  the  analytic  solution.  The  error  depends  on  the  grid  resolution  mod  3  since 
this  value  determines  how  the  structure  aligns  the  domain. 


to  resolution.  For  higher  resolutions,  the  fluid  time  step  restriction  will  dominate.  We  note  that  since  all  of 
our  examples  converge  to  nonzero  velocities  under  refinement,  we  effectively  refine  with  At  oc  Ax.  Thus,  the 
first  order  convergence  demonstrated  in  our  examples  implies  that  our  scheme  is  at  least  first  order  accurate 
in  space  and  time.  We  do  not  test  convergence  order  separately  for  space  and  time,  since  we  do  not  expect 
to  achieve  better  than  first  order. 

9.1.  Analytic  example 

Figure  3  shows  a  simple  setup  designed  to  mimic  an  infinite  channel  with  a  rigid  structure  flowing  through 
it  and  centered  in  the  channel.  The  density  of  the  fluid  is  p  =  100 kgm~2,  and  the  mass  of  the  rigid  body  is 
M  —  150 kg.  The  height  of  the  fluid  domain  is  h  =  1  m,  and  the  total  domain  width  of  1  m  is  divided  into 
two  fluid  region  with  width  w  =  ^ m  with  the  central  occupied  by  the  structure.  The  fluid  viscosity 
is  p  =  lOOfcgs-1.  Both  the  fluid  and  structure  experience  gravity  g  =  —9.8 ms-1.  We  elongate  the  rigid 
body  to  extend  slightly  beyond  the  domain  and  periodically  reset  its  position,  providing  the  illusion  of  an 
infinitely  long  rigid  body  and  allowing  the  simulation  to  be  run  long  enough  to  reach  steady  state.  Altering 
the  dimensions  of  the  structure  does  not  affect  the  analytic  solution,  provided  its  mass  is  not  altered. 

We  derive  the  analytic  solution  by  considering  the  structure  to  be  a  boundary  condition  to  the  fluid, 
and  then  solving  for  the  structure’s  velocity  using  conservation  of  energy.  The  analytic  velocity  profile  for 
the  left  fluid  region  f l  =  [0,  u>]  x  [0,  h]  is  it  =  0  and  v  =  §^x(x  —  w)  +  ^x,  where  the  v  is  the  vertical  fluid 
velocity  component  and  vs  is  the  vertical  component  of  the  structure  velocity.  The  fluid  experiences  no 
pressure  gradient.  The  only  unknown  remaining  is  vs,  the  velocity  of  the  structure.  Taking  into  account  the 
symmetry  of  the  system,  conservation  of  energy  implies  that  the  viscous  dissipation  equals  the  total  work 
applied  by  gravity,  so  that  2  fn(r  ■  V)u  dV  =  2  fn  pgv  dV  +  Mgvs.  Solving  this  equation  yields  the  analytic 
solution  vs  =  —  ( M  +  pwh )  . 

9.2.  Rigid  cylinder  falling  under  gravity 

Motivated  by  [30,  1],  a  rigid  cylinder  is  allowed  to  fall  under  an  externally  applied  gravitational  force. 
The  fluid  domain  has  dimensions  2 L  x  8 L1  and  the  cylinder  has  radius  r.  If  the  fluid  density  is  pp  and  the 
solid  density  is  ps,  then  the  terminal  velocity  is  [30] 

”'™  =  -  Hr1  (- ln  (z)  -  °-9157 + L7244(i Y  -  L7302(i  Y) ■ 

It  is  instructive  to  note  the  origin  of  this  analytic  solution.  Consider  a  fluid  domain  consisting  of  an 
infinite  2D  channel  with  parallel  walls  separated  by  a  distance  of  2 L.  In  the  middle  of  this  channel  is  a 
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(a)  Navier- Stokes,  g  —  0-1  kgs  1 


(c)  Stokes  flow,  g  =  0.1  kgs  1 


Figure  4:  Steady  state  velocity  field. 


Time  (s) 


Time  (s) 


(a)  Navier- Stokes,  g  =  0.1  kgs  1 


(b)  Navier- Stokes,  g  =  0.2 kgs  1 


(c)  Navier- Stokes,  g,  =  0.5 kgs 


(d)  Stokes  flow,  g  =  0.1  kgs 


Figure  5:  Velocity  of  a  falling  cylinder  in  a  channel  over  time. 
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Equations 

Viscosity 

(kgs-1) 

Analytic 

Terminal  Velocities  (ms  *) 
40  x  160  80  x  320 

160  x  640 

Convergence  order 

N-S 

0.1 

-0.3501 

-0.1716 

-0.1882 

-0.1966 

- 

- 

N-S 

0.2 

-0.1750 

-0.1254 

-0.1359 

-0.1417 

- 

- 

N-S 

0.5 

-0.07002 

-0.05983 

-0.06438 

-0.06721 

0.85 

1.01 

N-S 

1.0 

-0.03501 

-0.03052 

-0.03264 

-0.03399 

0.92 

1.22 

N-S 

2.0 

-0.01750 

-0.01533 

-0.01635 

-0.01702 

0.91 

1.26 

N-S 

5.0 

-0.007002 

-0.006148 

-0.006547 

-0.006828 

0.91 

1.38 

N-S 

10.0 

-0.003501 

-0.003013 

-0.003279 

-0.003417 

1.14 

1.40 

Stokes 

0.1 

-0.3501 

-0.2945 

-0.3241 

-0.3390 

1.10 

1.23 

Figure  6:  Terminal  velocity  of  a  falling  cylinder  and  convergence  orders. 


cylinder  with  radius  a  located  in  the  middle  of  the  channel  and  moving  with  a  constant  velocity  vs  along  the 
axis  of  the  channel.  Under  the  simplifying  assumption  of  Stokes  flow,  the  resistive  force  experienced  by  the 
cylinder  is  approximately  [15] 

F  _  _ 4:TTfXVg _ 

-  In  (f)  -0.9157+ 1.7244(f)2  -  1.7302(f)4' 

Note  that  because  of  the  Stokes  flow  assumption,  fluid  advection  and  movement  of  the  cylinder  are  assumed 
to  be  insignificant  and  dropped.  The  steady  state  solid  velocity  is  obtained  by  balancing  the  resistive  force 
against  the  net  gravitational  force  F  =  (ps  —  pf)gnr2  on  the  solid  due  to  gravity  (taking  into  account 
buoyancy) . 

We  use  the  same  parameters  as  [30]:  g  =  9.8 ?ns_1,  r  =  5  x  10-3  rn,  L  =  2  x  10~2?n,  p  =  0.1  kgs-1, 
Pf  =  lx  103kgm~ 2,  and  ps  =  2  x  10  3kgm~2.  The  top  of  the  fluid  domain  has  a  zero  pressure  Dirichlet 
boundary  condition,  and  the  sides  and  bottom  are  no-slip  walls.  Our  results  are  shown  in  Figure  9.2.  We 
note  that  our  results  do  not  agree  with  those  in  [30],  and  we  note  that  their  velocity  profile  does  appear  to 
reach  a  terminal  velocity.  We  converge  to  a  value  significantly  lower  than  the  analytic  solution.  Noting  that 
the  analytic  solution  was  derived  under  the  assumption  of  Stokes  flow,  we  also  ran  with  higher  viscosities: 
p  =  0.2  kgs-1  (Figure  9.2),  p,  =  0.5  kc /s~1  (Figure  9.2),  p  =  1.0  kgs-1,  p  =  2.0  kg  s^1 ,  p  =  5.0  kgs-1, 
and  p  =  10.0  kgs-1.  We  obtain  close  agreement  with  the  analytic  solution  for  p  =  0.5  kg  s^1  and  higher 
viscosities.  We  also  repeated  the  the  original  test  p  =  0.1  kgs-1  using  Stokes  flow  and  converge  to  the 
analytic  terminal  velocity  as  shown  in  Figure  9.2.  The  results  from  all  of  these  simulations  are  summarized 
in  Figure  6,  and  the  steady  state  velocities  are  shown  in  Figure  4. 

The  velocity  profile  shown  in  Figure  7  of  [1]  does  show  a  faster  initial  increase  in  velocity  followed  by  the 
velocity  profile  reaching  a  plateau,  which  is  consistent  with  the  behavior  we  observe.  Their  method,  like  ours, 
is  first  order  accurate,  and  their  figure  shows  two  velocity  profiles  for  each  of  two  methods  corresponding 
to  a  factor-of-four  refinement.  This  suggests  that  the  close  fit  shown  in  their  highest  resolution  may  be 
coincidental  and  that  their  scheme  will  actually  converge  to  a  value  somewhat  higher  than  the  predicted 
value. 

9.3.  Laminar  flow  past  a  cylinder 

In  this  example,  we  consider  a  simple  one-way  coupled  setup  and  directly  compare  compare  our  technique 
against  a  different  scheme.  For  the  comparison  scheme,  we  use  standard  Chorin  splitting  [8].  We  then  choose 
a  sample  point  at  a  dynamic  location  and  time  and  demonstrate  that  the  two  schemes  converge  to  the  same 
velocity  at  that  point.  This  demonstrates  that  in  the  simpler  case  where  structures  are  only  one-way  coupled 
to  the  fluid,  we  are  consistent  with  an  existing  scheme  from  the  literature. 

For  this  comparison  we  use  laminar  flow  across  a  fixed  cylinder.  The  domain  is  [0m,4m]  x  [0  rn.  4  rn] 
with  a  cylinder  of  radius  0.5  m  located  at  (2,  2).  The  cylinder  has  no-slip  boundary  conditions.  The  top  and 
bottom  walls  have  slip  boundary  conditions.  The  right  wall  has  the  condition  p  =  0,  and  the  left  wall  has 
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(0,0) 


(4,4) 


(a)  Domain  configuration  (units  in  m) 


(b)  Velocity  field 


loglO(resolution) 

(c)  Convergence  of  ui  at  A. 


(d)  Convergence  of  112  at  A. 


Figure  7:  Laminar  flow  around  a  fixed  cylinder. 


19 


the  boundary  condition  Ui  =  1  and  a  slipping  boundary  condition  for  U2.  The  density  is  p  =  1  kgm~2,  and 
the  viscosity  is  /i  =  0.1  kg  s~1.  The  initial  fluid  velocity  is  u  =  0.  The  setup  is  illustrated  in  Figure  7(a). 
The  fluid  velocity  is  sampled  at  the  point  (3,  2.25)  at  time  0.5s  using  linear  interpolation.  The  fluid  velocity 
is  far  from  steady  state  at  this  time,  so  this  comparison  also  tests  dynamic  behavior. 

The  Chorin  splitting  implementation  at  resolution  640  x  640  is  taken  to  be  the  exact  solution.  The 
example  was  then  run  using  the  proposed  method  at  resolutions  40  x  40,  80  x  80,  160  x  160,  and  320  x  320. 
The  differences  for  each  velocity  component  are  plotted  in  Figure  7(c)  and  Figure  7(d).  The  convergence 
order  for  Ui  is  0.89,  and  the  convergence  order  for  U2  is  1.07.  The  final  velocity  profile  is  shown  in  Figure  7(b). 

9-4-  Vortex  shedding 

In  this  example  we  demonstrate  Karman  vortex  shedding  around  a  cylinder  and  compare  the  dynamic 
behavior  of  our  method  in  the  case  of  low  viscosity  and  one-way  coupling  against  the  results  that  others  have 
obtained  and  against  our  implementation  of  Chorin  splitting.  Our  setup  is  based  on  that  of  [17,  6],  except 


(a)  Proposed  method 


(b)  Chorin  splitting 

Figure  8:  Vorticity  contours  demonstrating  vortex  shedding.  Contours  are  provided  for  the  range  [—2.5,  2.5]  in  increments  of 
0.25. 
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loglO(resolution) 

(b)  Position  error  ||x  —  xex 


Figure  9:  Velocity  profile  and  position  error  for  the  top  right  particle  of  the  deforming  beam  at  time  1  s.  The  position  of  the 
particle  at  the  highest  resolution  is  taken  to  be  the  exact  value  xex. 


that  we  employ  slightly  different  boundary  conditions.  Our  domain  is  [—2.5  m,  15.0  m]  x  [— 3.5  m,  3.5  m]  with 
a  fixed  cylinder  of  radius  0.5  m  centered  at  the  origin.  The  top  and  bottom  of  the  domain  walls  have  slip 
boundary  conditions.  The  left  boundary  has  a  source  term  Ui  =  1ms-1  with  a  slipping  boundary  condition 
for  U2.  The  right  boundary  has  the  boundary  condition  p  =  0.  The  flow  has  Re  =  100  with  p  =  1  kgm-2  and 
p  =  0.01  kgs-1.  We  used  a  400  x  160  grid  for  our  proposed  method  and  for  Chorin  splitting.  Our  vorticity 
contours  for  our  proposed  method  and  Chorin  splitting  are  shown  in  Figure  8.  Our  vorticity  contours  exhibit 
the  same  characteristics  as  those  shown  in  [17,  6]  as  well  as  being  quite  similar  to  each  other.  We  note  that 
the  contour  plot  for  the  proposed  has  an  additional  zero  contour  near  the  left  side.  The  vorticity  in  this  region 
is  quite  small,  reaching  its  peak  at  roughly  1%  of  the  separation  between  two  contour  lines.  At  resolution 
200  x  80,  the  vorticity  peaks  at  roughly  3%  of  the  separation  between  two  contour  lines,  indicating  that  the 
extra  contour  line  will  go  away  under  refinement. 

9. 5.  Flow  past  deformable  beam 

This  problem  deals  with  a  deformable  beam  immersed  in  fluid,  similar  to  problems  investigated  by 
[3,  31,  29,  30,  9,  27,  19,  32],  The  fluid  domain  was  [0m,2ro]  x  [0m,  1  m]  with  velocity  boundary  conditions 
of  0.25  ms-1  on  the  left  side,  slip  walls  on  the  top  and  bottom,  and  a  p  =  0  boundary  condition  on  the  right 
wall.  The  beam  has  no-slip  boundary  conditions.  The  fluid  density  is  100  kgm-2,  and  the  coefficient  of 
dynamic  viscosity  is  p  =  1  kgs-1.  The  structure  is  discretized  as  a  triangle  mesh  using  a  mass-spring  model 
with  20  particles  in  the  vertical  direction  and  5  in  the  horizontal,  each  of  mass  0.1  kg  except  for  the  bottom 
row  of  particles,  which  are  treated  as  infinitely  massive.  The  mesh  is  initially  rectangular,  with  extent  from 
[0.95 m,  1.05 m]  x  [0m,0.5m].  Edge  springs  with  Young’s  modulus  E  =  166 kgm-1s-2  and  overdamping 
fraction  2  along  with  altitude  springs  with  Young’s  modulus  E  =  400  kgm-1s-2  and  overdamping  fraction 
2  were  used  as  the  structure  constitutive  model.  (If  considered  in  isolation,  a  spring  can  be  made  critically 
damped  by  choosing  a  specific  damping  constant.  When  we  refer  to  an  overdamping  fraction  of  2,  we  mean 
that  the  damping  coefficient  on  the  spring  is  twice  the  value  that  would  result  in  a  critically  damped  spring. 
Isolated  springs  with  an  overdamping  fraction  greater  than  one  are  overdamped,  and  isolated  springs  with  an 
overdamping  fraction  less  than  one  are  underdamped.  The  analogy  does  not  extend  to  systems  of  springs, 
but  we  find  this  approach  to  be  more  intuitive  for  describing  the  amount  of  damping  present.)  Altitude 
springs  [21]  create  a  spring  between  each  particle  of  a  triangle  and  a  virtual  node  projected  onto  the  opposite 
face  and  are  used  to  give  volumetric  structure  and  prevent  inversion. 

Note  for  the  purposes  of  our  convergence  analysis  we  take  the  structure  to  be  a  discrete  system  and 
refine  only  the  fluid  grid  resolution.  We  examine  convergence  with  fluid  grid  resolution  and  time  step  size  in 
Figure  9(b).  The  error  is  computed  assuming  the  simulation  with  resolution  256  x  128  is  the  exact  solution. 
We  plot  the  convergence  of  ||x  —  xea;||,  where  x.ex  is  the  position  in  the  exact  solution  configuration,  for  the 
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(c)  t  =  0.500  s 


(d)  t  =  1.00  s 


Figure  10:  Streamlines  in  the  oscillating  circle  example.  The  fluid  grid  resolution  is  128  X  128  and  the  structure  is  composed  of 
2400  triangles. 


top  right  particle  of  the  beam  at  time  1  s  in  Figure  9(b).  The  convergence  order  is  1.42.  The  velocity  profile 
is  shown  in  Figure  9(a). 

9.6.  Oscillating  circle 

We  have  implemented  the  oscillating  disk  example  from  [32] ,  where  a  circular  deformable  body  is  placed 
within  a  fluid  domain  of  dimensions  [0  to,  1  to]  x  [0  to,  1  to]  with  periodic  boundary  conditions  and  an  initial 
velocity  defined  for  both  the  fluid  and  the  structure  by  the  stream  function  =  iposin{kxx)sin{kyy ),  with 
if) o  =  0.05to2s-1  and  kx  =  ky  =  2'irm~1.  We  used  a  neo-Hookean  constitutive  model  as  described  in  [5]  for 
the  structure  with  Poisson’s  ratio  v  =  0.475,  shear  modulus  G  =  1  Pa,  and  Rayleigh  damping  coefficient 
0.005.  The  structure  and  fluid  densities  were  set  to  p  =  1  kgm~2  and  the  coefficient  of  dynamic  viscosity 
was  y  =  0.001  kgs-1.  We  ran  the  simulation  for  a  period  of  Is.  For  comparison  to  [32]  we  rasterize  the 
structure  velocities  onto  the  fluid  grid  and  display  the  streamlines  and  structure  deformation  (Figure  10). 
The  frequency  of  the  oscillator  is  similar  to  that  in  [32]  as  can  be  seen  from  the  energy  plots  in  Figure  11.  We 
note  that  our  simulation  is  more  damped  than  theirs.  Our  extra  damping  is  due  to  our  use  of  semi-Lagrangian 
advection  on  the  fluid  and  Rayleigh  damping  on  our  structure. 

10.  Conclusions  and  Future  Work 

We  have  presented  a  method  for  formulating  strongly  coupled  fluid  structure  interaction  problems  as  a 
symmetric  positive  definite  linear  system.  The  method  is  based  on  factoring  the  structure  damping  matrix 
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Figure  11:  Energy  over  time  for  the  oscillating  circle  shown  in  Figure  10. 


and  performing  a  change  of  variables  in  the  structure  unknowns.  It  can  be  applied  to  arbitrary  linearized 
structure  constitutive  models,  but  requires  a  semi-implicit  time  discretization  for  the  structure  equations 
where  elastic  terms  are  handled  explicitly  and  damping  terms  are  handled  implicitly.  For  the  special  case 
of  only  rigid  structures  and  no  other  rigid-rigid  forces  we  recover  exactly  the  method  of  [4] .  The  numerical 
experiments  on  problems  common  in  the  fluid  structure  interaction  literature  demonstrate  that  the  method 
is  convergent  with  grid  resolution. 

Second  order  accuracy  for  fluid  structure  interaction  has  been  achieved  by  some  authors,  such  as  the 
formally  second  order  methods  [13,  12].  We  note,  however,  that  second  order  accuracy  came  well  after  the 
original  method  and  that  the  additional  accuracy  comes  at  the  cost  of  additional  complexity.  Instead  of 
attempting  to  address  second  order  accuracy  in  this  paper,  we  leave  the  problem  of  higher  order  accuracy 
for  future  work. 

Compared  to  standard  Chorin  splitting  (without  structures),  our  evolution  will  be  significantly  slower. 
Chorin  splitting  with  implicit  viscosity  requires  the  solution  of  one  well-conditioned  viscosity  system  per 
dimension  and  a  single  poorly- conditioned  Poisson  system  for  which  preconditioning  has  been  well  studied. 
We  replace  these  with  a  single  linear  system  that,  due  to  the  Poisson  matrix  it  contains,  is  also  poorly 
conditioned.  We  have  observed  its  conditioning  to  be  significantly  worse  at  high  viscosity  than  low  viscosity, 
and  we  have  also  observed  poor  conditioning  to  result  from  underresolved  regions.  In  particular,  we  have 
observed  poor  convergence  when  the  structure  is  much  coarser  than  the  fluid  grid  and  we  add  extra  coupling 
faces  to  W  as  described  in  Section  7,  possibly  due  to  overconstraint.  Increasing  the  resolution  of  the 
structure  improves  convergence  in  these  cases,  but  a  more  careful  treatment  of  boundary  conditions  might 
avoid  the  issue  entirely.  Conditioning  with  deformable  structures  becomes  worse  with  lower  density  and 
higher  damping.  The  conditioning  improves  significantly  when  all  structures  are  one-way  coupled. 

The  system  is  symmetric  positive  definite,  making  it  suitable  for  solution  with  methods  like  conjugate 
gradient  and  generally  easier  to  precondition  and  solve  than  to  other  types  of  linear  systems.  The  solution 
of  this  system  is  still  expensive  in  our  implementation  as  we  have  not  developed  an  effective  preconditioner 
for  it  and  have  not  made  significant  effort  to  optimize  our  routines.  We  leave  the  problem  of  finding  a  good 
preconditioner  for  future  work.  Solving  the  linear  systems  dominates  the  cost  of  our  time  evolution.  The 
algorithm  compares  favorably  with  [25],  which  has  an  indefinite  system  that  we  have  observed  to  converge 
much  slower  in  practice. 
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